rm(list=ls())
## Set working directory
setwd("/Users/afanlo/Documents/Projects/Kroenig Rep")
library(dplyr)
library(MASS)
library(xtable)
library(stargazer)
library(lmtest)
library(sandwich)

## Upload Kroenig Data
cutoff_df<-read.csv("DF_WithCutoffs.csv")
cutoff_df<-cutoff_df[,c(3:9, 29:148)]
dataframe<-read.csv("data_may2020.csv")
df<-merge(dataframe, cutoff_df, by=c("ccode", "ccode2", "crisno", "year", "crisname",
                                     "abbrev1", "abbrev2"))

##Changing back to ICB coding for Korean War
df$victory<-ifelse(df$crisno==133 & df$ccode==365,  1, df$victory)

## Changing back to ICB Coding of 1969 Sino-Soviet Border War
df$victory<-ifelse(df$crisno==231, 0, df$victory)

#Dropping France Berlin Wall Dyad
df$frn_berlin<-ifelse(df$crisno==185 & df$ccode==220 | df$crisno==185 & df$ccode2==220, 1, 0)
df<- subset(df, df$frn_berlin==0)

df$inferiority<-ifelse(df$superiority==0, 1, 0)

### Creating a function to compute clustered standard errors
vcovCluster <- function(
    model,
    cluster
)
{
  require(sandwich)
  require(lmtest)
  
  cluster = as.factor(cluster)
  
  if (!is.null(model$na.action)){
    omit.rows = model$na.action
    cluster = cluster[-omit.rows]
  }
  if (nrow(model.matrix(model)) != length(cluster)){
    stop("something's not working: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap or randomization inference instead).")
  }
  
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- na.omit(apply(estfun(model), 2, function(x) tapply(x, cluster, sum)))
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

## Models Begin


## 27 1s
logit_1<-glm(data=df, victory~superiority+proximity  + polity21 + 
               capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test1<-coeftest(logit_1, vcov=vcovCluster(model=logit_1, cluster=factor(df$crisisdyad)))


logit_1.15<-glm(data=df, victory~onep15cut+op15_par+proximity  + polity21 + 
                  capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test1.15<-coeftest(logit_1.15, vcov=vcovCluster(model=logit_1.15, cluster=factor(df$crisisdyad)))

### 26 1s
logit_1.2<-glm(data=df, victory~onep2cut+onep2_par+ proximity  + polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test1.2<-coeftest(logit_1.2, vcov=vcovCluster(model=logit_1.2, cluster=factor(df$crisisdyad)))

#### 25 1s
table(df$onep45cut)

logit_1.45<-glm(data=df, victory~onep45cut+op45_par+proximity +  polity21 + 
                  capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test1.45<-coeftest(logit_1.45, vcov=vcovCluster(model=logit_1.45, cluster=factor(df$crisisdyad)))

## 24 1s
logit_1.55<-glm(data=df, victory~onep55cut+op55_par+proximity +  polity21 + 
                  capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test1.55<-coeftest(logit_1.55, vcov=vcovCluster(model=logit_1.55, cluster=factor(df$crisisdyad)))


## 23 1s
logit_1.65<-glm(data=df, victory~onep65cut+op65_par+proximity +  polity21 + 
                  capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test1.65<-coeftest(logit_1.65, vcov=vcovCluster(model=logit_1.65, cluster=factor(df$crisisdyad)))

## 22 1s
logit_1.85<-glm(data=df, victory~onep85cut+op85_par+proximity +  polity21 + 
                  capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test1.85<-coeftest(logit_1.85, vcov=vcovCluster(model=logit_1.85, cluster=factor(df$crisisdyad)))


## 21 1s
logit_2<-glm(data=df, victory~twocut+two_par+proximity +  polity21 + 
               capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test2<-coeftest(logit_2, vcov=vcovCluster(model=logit_2, cluster=factor(df$crisisdyad)))


### 20 1s
logit_2.5<-glm(data=df, victory~twop5cut+tp5_par+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test2.5<-coeftest(logit_2.5, vcov=vcovCluster(model=logit_2.5, cluster=factor(df$crisisdyad)))


### 19 1s
logit_3<-glm(data=df, victory~threecut+three_par+proximity +  polity21 + 
               capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test3<-coeftest(logit_3, vcov=vcovCluster(model=logit_3, cluster=factor(df$crisisdyad)))


### 18 1s
logit_3.4<-glm(data=df, victory~threep4cut+threep4_par+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test3.4<-coeftest(logit_3.4, vcov=vcovCluster(model=logit_3.4, cluster=factor(df$crisisdyad)))

### 17 1s
logit_3.5<-glm(data=df, victory~threep5cut+threep5_par+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test3.5<-coeftest(logit_3.5, vcov=vcovCluster(model=logit_3.5, cluster=factor(df$crisisdyad)))


### 16 1s 
logit_4<-glm(data=df, victory~fourcut+four_par+proximity +  polity21 + 
               capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test4<-coeftest(logit_4, vcov=vcovCluster(model=logit_4, cluster=factor(df$crisisdyad)))

### 15 1s
logit_7<-glm(data=df, victory~sevencut+seven_par+proximity +  polity21 + 
               capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test7<-coeftest(logit_7, vcov=vcovCluster(model=logit_7, cluster=factor(df$crisisdyad)))

### 14 1s
logit_9<-glm(data=df, victory~ninecut+nine_par+proximity +  polity21 + 
               capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test9<-coeftest(logit_9, vcov=vcovCluster(model=logit_9, cluster=factor(df$crisisdyad)))

### 13 1s
logit_10<-glm(data=df, victory~tencut+ten_par+proximity +  polity21 + 
                capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test10<-coeftest(logit_10, vcov=vcovCluster(model=logit_10, cluster=factor(df$crisisdyad)))

## 12 1s

logit_11<-glm(data=df, victory~elevencut+elevencut_par+proximity +  polity21 + 
                capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test11<-coeftest(logit_11, vcov=vcovCluster(model=logit_11, cluster=factor(df$crisisdyad)))


### 11 1s
logit_20<-glm(data=df, victory~twentycut+twenty_par+proximity +  polity21 + 
                capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test20<-coeftest(logit_20, vcov=vcovCluster(model=logit_20, cluster=factor(df$crisisdyad)))

### 10 1s
logit_28<-glm(data=df, victory~twentyeightcut+twentyeight_par+proximity +  polity21 + 
                capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test28<-coeftest(logit_28, vcov=vcovCluster(model=logit_28, cluster=factor(df$crisisdyad)))


### 9 1s
logit_30<-glm(data=df, victory~thirtycut+thirty_par+proximity +  polity21 + 
                capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test30<-coeftest(logit_30, vcov=vcovCluster(model=logit_30, cluster=factor(df$crisisdyad)))

### 8 1s
logit_40<-glm(data=df, victory~fortycut+forty_par+proximity +  polity21 + 
                capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test40<-coeftest(logit_40, vcov=vcovCluster(model=logit_40, cluster=factor(df$crisisdyad)))

### 7 1s
logit_50<-glm(data=df, victory~fiftycut+fifty_par+proximity +  polity21 + 
                capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test50<-coeftest(logit_50, vcov=vcovCluster(model=logit_50, cluster=factor(df$crisisdyad)))

### 6 1s
logit_100<-glm(data=df, victory~hundredcut+hundred_par+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test100<-coeftest(logit_100, vcov=vcovCluster(model=logit_100, cluster=factor(df$crisisdyad)))

### 5 1s 
logit_500<-glm(data=df, victory~fivhuncut+fivhun_par+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test500<-coeftest(logit_500, vcov=vcovCluster(model=logit_500, cluster=factor(df$crisisdyad)))

### 4 1s 

logit_700<-glm(data=df, victory~sevenhuncut+sevenhuncut_par+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test700<-coeftest(logit_700, vcov=vcovCluster(model=logit_700, cluster=factor(df$crisisdyad)))

### 3 1s 
logit_1000<-glm(data=df, victory~thoucut+thou_par+proximity +  polity21 + 
                  capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test1000<-coeftest(logit_1000, vcov=vcovCluster(model=logit_1000, cluster=factor(df$crisisdyad)))


# 2 1s
logit_1300<-glm(data=df, victory~thirhuncut+thirhun_par+proximity +  polity21 + 
                  capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
test1300<-coeftest(logit_1300, vcov=vcovCluster(model=logit_1300, cluster=factor(df$crisisdyad)))

#1 1s
table(df$fifthuncut)
##


## Preparing Data Frames for First Differences

superior<-as.data.frame(t(as.matrix(c(1, mean(df$proximity), mean(df$polity21),
                                      mean(df$capshare), mean(df$strikea), mean(df$tpop_1), 
                                      mean(df$viol), mean(df$cris_ave)))))
colnames(superior)<-c("superiority", "proximity",
                      "polity21", "capshare", "strikea", "tpop_1",
                      "viol", "cris_ave")

inferior<-as.data.frame(t(as.matrix(c(0, mean(df$proximity), mean(df$polity21),
                                      mean(df$capshare), mean(df$strikea), mean(df$tpop_1), 
                                      mean(df$viol), mean(df$cris_ave)))))
colnames(inferior)<-c("superiority", "proximity",
                      "polity21", "capshare", "strikea", "tpop_1",
                      "viol", "cris_ave")


######

superior1<-as.data.frame(t(as.matrix(c(1, 0, mean(df$proximity), mean(df$polity21),
                                       mean(df$capshare), mean(df$strikea), mean(df$tpop_1), 
                                       mean(df$viol), mean(df$cris_ave)))))
colnames(superior1)<-c("superiority", "parity", "proximity",
                       "polity21", "capshare", "strikea", "tpop_1",
                       "viol", "cris_ave")

inferior1<-as.data.frame(t(as.matrix(c(0, 0, mean(df$proximity),mean(df$polity21),
                                       mean(df$capshare), mean(df$strikea), mean(df$tpop_1), 
                                       mean(df$viol), mean(df$cris_ave)))))
colnames(inferior1)<-c("superiority", "parity", "proximity",
                       "polity21", "capshare", "strikea", "tpop_1",
                       "viol", "cris_ave")

first_diffs<-matrix(NA, 27, 8)
first_diffs[,1]<-c("logit_1", "logit_1.15", "logit_1.2",
                   "logit_1.45", "logit_1.55", "logit_1.65", "logit_1.85",
                   "logit_2", "logit_2.5", "logit_3", "logit_3.4", "logit_3.5",
                   "logit_4", "logit_7", 
                   "logit_9", "logit_10", "logit_11", "logit_20", "logit_28", "logit_30", 
                   "logit_40",   "logit_50",  
                   "logit_100", "logit_500", "logit_700", "logit_1000", "logit_1300")
first_diffs[,2]<-seq(from=28, to=2, by=-1)
first_diffs[,3]<-c(names(logit_1$coefficients[2]), names(logit_1.15$coefficients[2]), names(logit_1.2$coefficients[2]),
                   names(logit_1.45$coefficients[2]), names(logit_1.55$coefficients[2]), names(logit_1.65$coefficients[2]), names(logit_1.85$coefficients[2]),
                   names(logit_2$coefficients[2]), names(logit_2.5$coefficients[2]), names(logit_3$coefficients[2]), names(logit_3.4$coefficients[2]), 
                   names(logit_3.5$coefficients[2]),
                   names(logit_4$coefficients[2]), names(logit_7$coefficients[2]), 
                   names(logit_9$coefficients[2]), names(logit_10$coefficients[2]), names(logit_11$coefficients[2]), names(logit_20$coefficients[2]), names(logit_28$coefficients[2]), 
                   names(logit_30$coefficients[2]), 
                   names(logit_40$coefficients[2]),   names(logit_50$coefficients[2]),  
                   names(logit_100$coefficients[2]), names(logit_500$coefficients[2]), names(logit_700$coefficients[2]), names(logit_1000$coefficients[2]), names(logit_1300$coefficients[2]))  

first_diffs[,4]<-c(NA, names(logit_1.15$coefficients[3]), names(logit_1.2$coefficients[3]),
                   names(logit_1.45$coefficients[3]), names(logit_1.55$coefficients[3]), names(logit_1.65$coefficients[3]), names(logit_1.85$coefficients[3]),
                   names(logit_2$coefficients[3]), names(logit_2.5$coefficients[3]), names(logit_3$coefficients[3]), names(logit_3.4$coefficients[3]), 
                   names(logit_3.5$coefficients[3]),
                   names(logit_4$coefficients[3]), names(logit_7$coefficients[3]), 
                   names(logit_9$coefficients[3]), names(logit_10$coefficients[3]), names(logit_11$coefficients[3]), names(logit_20$coefficients[3]), names(logit_28$coefficients[3]), 
                   names(logit_30$coefficients[3]), 
                   names(logit_40$coefficients[3]),   names(logit_50$coefficients[3]),  
                   names(logit_100$coefficients[3]), names(logit_500$coefficients[3]), names(logit_700$coefficients[3]), names(logit_1000$coefficients[3]), names(logit_1300$coefficients[3]))  

mods<-list(logit_1, logit_1.15, logit_1.2,
           logit_1.45, logit_1.55, logit_1.65, logit_1.85,
           logit_2, logit_2.5, logit_3, logit_3.4, logit_3.5,
           logit_4, logit_7, 
           logit_9, logit_10, logit_11, logit_20, logit_28, logit_30, 
           logit_40,   logit_50,  
           logit_100, logit_500, logit_700, logit_1000, logit_1300)
for (i in 2:27){
  sup_new<-superior1
  colnames(sup_new)<-c(first_diffs[i,3], first_diffs[i, 4], colnames(superior1)[3:9]) 
  inf_new<-inferior1
  colnames(inf_new)<-c(first_diffs[i,3], first_diffs[i, 4], colnames(inferior1)[3:9]) 
  mod_new<-mods[[i]]
  first_diffs[i,5]<-predict(mod_new, sup_new, type="response")
  first_diffs[i,6]<-predict(mod_new, inf_new, type="response")
  sup_new<-NA
  inf_new<-NA
}

first_diffs[1,5]<-predict(logit_1, superior, type="response")
first_diffs[1,6]<-predict(logit_1, inferior, type="response")

first_diffs[,7]<-as.numeric(first_diffs[,5])-as.numeric(first_diffs[,6])

## Determining Which Coefs on Sup are Significant, According to Cluster-Robust SEs
sig<-matrix(NA, 27, 4)
sig[,1]<-c("logit_1", "logit_1.15", "logit_1.2",
           "logit_1.45", "logit_1.55", "logit_1.65", "logit_1.85",
           "logit_2", "logit_2.5", "logit_3", "logit_3.4", "logit_3.5",
           "logit_4", "logit_7", 
           "logit_9", "logit_10", "logit_11", "logit_20", "logit_28", "logit_30", 
           "logit_40",   "logit_50",  
           "logit_100", "logit_500", "logit_700", "logit_1000", "logit_1300")
sig[,2]<-c(test1[2,4], test1.15[2,4], test1.2[2,4],
           test1.45[2,4], test1.55[2,4], test1.65[2,4], test1.85[2,4],
           test2[2,4], test2.5[2,4], test3[2,4], test3.4[2,4], test3.5[2,4],
           test4[2,4], test7[2,4], 
           test9[2,4], test10[2,4], test11[2,4], test20[2,4], test28[2,4], test30[2,4], 
           test40[2,4],   test50[2,4],  
           test100[2,4], test500[2,4], test700[2,4], test1000[2,4], test1300[2,4])
sig<-as.data.frame(sig)
sig[,2]<-as.numeric(sig[,2])
sig[,3]<-ifelse(sig[,2]<0.1, 1, 0)
sig[,4]<-seq(1, 27, 1)

#### Figure 4 from Appendix C- With superiority using cluster-robust standard errors 
dev.off()
par(mar=c(8.1, 4.1, 4.1, 4.1), xpd=T)
plot(first_diffs[,7], ylim=c(-1,1), xaxt="n", xlab=c("Nuclear Ratio: Asymmetric Superiority Threshold"), 
     ylab=c("First Difference"), main=c("Difference in Prob. of Victory: Superior vs. Inferior"))
axis(1, at=1:27, labels=c("1", "1.15", "1.2",
                          "1.45", "1.55", "1.65", "1.85",
                          "2", "2.5", "3", "3.4", "3.5",
                          "4", "7", 
                          "9", "10", "11", "20", "28", "30", 
                          "40",   "50",  
                          "100", "500", "700", "1000", "1300"), las=2)
legend("bottomright", inset=c(0, -0.25),legend=c("Significant Effect, Cluster-Robust"), pch=19, 
       cex=0.5, pt.cex=1)
abline(h=0, xpd=F)

first_diffs[,8]<-seq(1,27)
points(first_diffs[1,8],first_diffs[1,7] ,  pch=19)
points(first_diffs[2,8],first_diffs[2,7] ,  pch=19)
points(first_diffs[4,8],first_diffs[4,7] ,  pch=19)
points(first_diffs[5,8],first_diffs[5,7] ,  pch=19)
points(first_diffs[6,8],first_diffs[6,7] ,  pch=19)
points(first_diffs[20,8],first_diffs[20,7] ,  pch=19)
points(first_diffs[22,8],first_diffs[22,7] ,  pch=19)
points(first_diffs[23,8],first_diffs[23,7] ,  pch=19)
points(first_diffs[24,8],first_diffs[24,7] ,  pch=19)
points(first_diffs[25,8],first_diffs[25,7] ,  pch=19)
points(first_diffs[26,8],first_diffs[26,7] ,  pch=19)
points(first_diffs[27,8],first_diffs[27,7] ,  pch=19)


########
## RANDOMIZATION INFERENCE 
######

rm(list=setdiff(ls(), "first_diffs"))
set.seed(02303)

library(dplyr)
library(MASS)
library(xtable)
library(stargazer)
library(lmtest)
library(sandwich)
setwd("/Users/afanlo/Documents/Projects/Kroenig Rep")

## Upload Kroenig Data
cutoff_df<-read.csv("DF_WithCutoffs.csv")
cutoff_df<-cutoff_df[,c(3:9, 29:148)]
dataframe<-read.csv("data_may2020.csv")
df<-merge(dataframe, cutoff_df, by=c("ccode", "ccode2", "crisno", "year", "crisname",
                                     "abbrev1", "abbrev2"))

##Changing back to ICB coding for Korean War
df$victory<-ifelse(df$crisno==133 & df$ccode==365,  1, df$victory)

## Changing back to ICB Coding of 1969 Sino-Soviet Border War
df$victory<-ifelse(df$crisno==231, 0, df$victory)

#Dropping France Berlin Wall Dyad
df$frn_berlin<-ifelse(df$crisno==185 & df$ccode==220 | df$crisno==185 & df$ccode2==220, 1, 0)
df<- subset(df, df$frn_berlin==0)


df$inferiority<-ifelse(df$superiority==0, 1, 0)

bin_df1p15_ss<-glm(data=df, victory~onep15cut+op15_par+proximity +  polity21 + 
                     capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
sims<-10000

coef.1p15.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$onep15cut)
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.1p15.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.1p15.ss<-cbind(coef.1p15.ss, fail)
coef.1p15.ss<-subset(coef.1p15.ss, coef.1p15.ss[,11]==0)

coef.1p151.ss<-matrix(NA, 1, 10)
coef.1p151.ss[1,]<-summary(bin_df1p15_ss)$coefficients[,3]

p.coef1p15.1<-mean(sapply(coef.1p15.ss[,1], FUN=function(x) x >= abs(coef.1p151.ss[1,1])))
p.coef1p15.2<-mean(sapply(coef.1p15.ss[,2], FUN=function(x) x >= abs(coef.1p151.ss[1,2])))
p.coef1p15.3<-mean(sapply(coef.1p15.ss[,3], FUN=function(x) x >= abs(coef.1p151.ss[1,3])))
p.coef1p15.4<-mean(sapply(coef.1p15.ss[,4], FUN=function(x) x >= abs(coef.1p151.ss[1,4])))
p.coef1p15.5<-mean(sapply(coef.1p15.ss[,5], FUN=function(x) x >= abs(coef.1p151.ss[1,5])))
p.coef1p15.6<-mean(sapply(coef.1p15.ss[,6], FUN=function(x) x >= abs(coef.1p151.ss[1,6])))
p.coef1p15.7<-mean(sapply(coef.1p15.ss[,7], FUN=function(x) x >= abs(coef.1p151.ss[1,7])))
p.coef1p15.8<-mean(sapply(coef.1p15.ss[,8], FUN=function(x) x >= abs(coef.1p151.ss[1,8])))
p.coef1p15.9<-mean(sapply(coef.1p15.ss[,9], FUN=function(x) x >= abs(coef.1p151.ss[1,9])))
p.coef1p15.10<-mean(sapply(coef.1p15.ss[,10], FUN=function(x) x >= abs(coef.1p151.ss[1,10])))

pvals_rand<-matrix(NA, 26, 10)
pvals_rand[1,]<-c(p.coef1p15.1, p.coef1p15.2, p.coef1p15.3, p.coef1p15.4, p.coef1p15.5,
                  p.coef1p15.6, p.coef1p15.7, p.coef1p15.8,  p.coef1p15.9, p.coef1p15.10)

discarded_iters<-matrix(NA, 27, 1)
discarded_iters[1,]<-10000-nrow(coef.1p15.ss)

### Cutoff 1.2
bin_df1p2_ss<-glm(data=df, victory~onep2cut+onep2_par+proximity +  polity21 + 
                    capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
sims<-10000

coef.1p2.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$onep2cut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.1p2.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.1p2.ss<-cbind(coef.1p2.ss, fail)
coef.1p2.ss<-subset(coef.1p2.ss, coef.1p2.ss[,11]==0)

coef.1p21.ss<-matrix(NA, 1, 10)
coef.1p21.ss[1,]<-summary(bin_df1p2_ss)$coefficients[,3]

p.coef1p2.1<-mean(sapply(coef.1p2.ss[,1], FUN=function(x) x >= abs(coef.1p21.ss[1,1])))
p.coef1p2.2<-mean(sapply(coef.1p2.ss[,2], FUN=function(x) x >= abs(coef.1p21.ss[1,2])))
p.coef1p2.3<-mean(sapply(coef.1p2.ss[,3], FUN=function(x) x >= abs(coef.1p21.ss[1,3])))
p.coef1p2.4<-mean(sapply(coef.1p2.ss[,4], FUN=function(x) x >= abs(coef.1p21.ss[1,4])))
p.coef1p2.5<-mean(sapply(coef.1p2.ss[,5], FUN=function(x) x >= abs(coef.1p21.ss[1,5])))
p.coef1p2.6<-mean(sapply(coef.1p2.ss[,6], FUN=function(x) x >= abs(coef.1p21.ss[1,6])))
p.coef1p2.7<-mean(sapply(coef.1p2.ss[,7], FUN=function(x) x >= abs(coef.1p21.ss[1,7])))
p.coef1p2.8<-mean(sapply(coef.1p2.ss[,8], FUN=function(x) x >= abs(coef.1p21.ss[1,8])))
p.coef1p2.9<-mean(sapply(coef.1p2.ss[,9], FUN=function(x) x >= abs(coef.1p21.ss[1,9])))
p.coef1p2.10<-mean(sapply(coef.1p2.ss[,10], FUN=function(x) x >= abs(coef.1p21.ss[1,10])))

pvals_rand[2,]<-c(p.coef1p2.1, p.coef1p2.2, p.coef1p2.3, p.coef1p2.4, p.coef1p2.5,
                  p.coef1p2.6, p.coef1p2.7, p.coef1p2.8,  p.coef1p2.9, p.coef1p2.10)

discarded_iters[2,]<-10000-nrow(coef.1p2.ss)

###Cutoff 1.45
bin_df1p45_ss<-glm(data=df, victory~onep45cut+op45_par+proximity +  polity21 + 
                     capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
sims<-10000

coef.1p45.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$onep45cut)
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.1p45.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.1p45.ss<-cbind(coef.1p45.ss, fail)
coef.1p45.ss<-subset(coef.1p45.ss, coef.1p45.ss[,11]==0)

coef.1p451.ss<-matrix(NA, 1, 10)
coef.1p451.ss[1,]<-summary(bin_df1p45_ss)$coefficients[,3]

p.coef1p45.1<-mean(sapply(coef.1p45.ss[,1], FUN=function(x) x >= abs(coef.1p451.ss[1,1])))
p.coef1p45.2<-mean(sapply(coef.1p45.ss[,2], FUN=function(x) x >= abs(coef.1p451.ss[1,2])))
p.coef1p45.3<-mean(sapply(coef.1p45.ss[,3], FUN=function(x) x >= abs(coef.1p451.ss[1,3])))
p.coef1p45.4<-mean(sapply(coef.1p45.ss[,4], FUN=function(x) x >= abs(coef.1p451.ss[1,4])))
p.coef1p45.5<-mean(sapply(coef.1p45.ss[,5], FUN=function(x) x >= abs(coef.1p451.ss[1,5])))
p.coef1p45.6<-mean(sapply(coef.1p45.ss[,6], FUN=function(x) x >= abs(coef.1p451.ss[1,6])))
p.coef1p45.7<-mean(sapply(coef.1p45.ss[,7], FUN=function(x) x >= abs(coef.1p451.ss[1,7])))
p.coef1p45.8<-mean(sapply(coef.1p45.ss[,8], FUN=function(x) x >= abs(coef.1p451.ss[1,8])))
p.coef1p45.9<-mean(sapply(coef.1p45.ss[,9], FUN=function(x) x >= abs(coef.1p451.ss[1,9])))
p.coef1p45.10<-mean(sapply(coef.1p45.ss[,10], FUN=function(x) x >= abs(coef.1p451.ss[1,10])))

pvals_rand[3,]<-c(p.coef1p45.1, p.coef1p45.2, p.coef1p45.3, p.coef1p45.4, p.coef1p45.5,
                  p.coef1p45.6, p.coef1p45.7, p.coef1p45.8,  p.coef1p45.9, p.coef1p45.10)

discarded_iters[3,]<-10000-nrow(coef.1p45.ss)


### Cutoff 1.55
bin_df1p55_ss<-glm(data=df, victory~onep55cut+op55_par+proximity +  polity21 + 
                     capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
sims<-10000

coef.1p55.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$onep55cut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.1p55.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.1p55.ss<-cbind(coef.1p55.ss, fail)
coef.1p55.ss<-subset(coef.1p55.ss, coef.1p55.ss[,11]==0)

coef.1p551.ss<-matrix(NA, 1, 10)
coef.1p551.ss[1,]<-summary(bin_df1p55_ss)$coefficients[,3]

p.coef1p55.1<-mean(sapply(coef.1p55.ss[,1], FUN=function(x) x >= abs(coef.1p551.ss[1,1])))
p.coef1p55.2<-mean(sapply(coef.1p55.ss[,2], FUN=function(x) x >= abs(coef.1p551.ss[1,2])))
p.coef1p55.3<-mean(sapply(coef.1p55.ss[,3], FUN=function(x) x >= abs(coef.1p551.ss[1,3])))
p.coef1p55.4<-mean(sapply(coef.1p55.ss[,4], FUN=function(x) x >= abs(coef.1p551.ss[1,4])))
p.coef1p55.5<-mean(sapply(coef.1p55.ss[,5], FUN=function(x) x >= abs(coef.1p551.ss[1,5])))
p.coef1p55.6<-mean(sapply(coef.1p55.ss[,6], FUN=function(x) x >= abs(coef.1p551.ss[1,6])))
p.coef1p55.7<-mean(sapply(coef.1p55.ss[,7], FUN=function(x) x >= abs(coef.1p551.ss[1,7])))
p.coef1p55.8<-mean(sapply(coef.1p55.ss[,8], FUN=function(x) x >= abs(coef.1p551.ss[1,8])))
p.coef1p55.9<-mean(sapply(coef.1p55.ss[,9], FUN=function(x) x >= abs(coef.1p551.ss[1,9])))
p.coef1p55.10<-mean(sapply(coef.1p55.ss[,10], FUN=function(x) x >= abs(coef.1p551.ss[1,10])))

pvals_rand[4,]<-c(p.coef1p55.1, p.coef1p55.2, p.coef1p55.3, p.coef1p55.4, p.coef1p55.5,
                  p.coef1p55.6, p.coef1p55.7, p.coef1p55.8,  p.coef1p55.9, p.coef1p55.10)

discarded_iters[4,]<-10000-nrow(coef.1p55.ss)



### Cutoff 1.65
bin_df1p65_ss<-glm(data=df, victory~onep65cut+op65_par+proximity +  polity21 + 
                     capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.1p65.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$onep65cut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.1p65.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.1p65.ss<-cbind(coef.1p65.ss, fail)
coef.1p65.ss<-subset(coef.1p65.ss, coef.1p65.ss[,11]==0)

coef.1p651.ss<-matrix(NA, 1, 10)
coef.1p651.ss[1,]<-summary(bin_df1p65_ss)$coefficients[,3]

p.coef1p65.1<-mean(sapply(coef.1p65.ss[,1], FUN=function(x) x >= abs(coef.1p651.ss[1,1])))
p.coef1p65.2<-mean(sapply(coef.1p65.ss[,2], FUN=function(x) x >= abs(coef.1p651.ss[1,2])))
p.coef1p65.3<-mean(sapply(coef.1p65.ss[,3], FUN=function(x) x >= abs(coef.1p651.ss[1,3])))
p.coef1p65.4<-mean(sapply(coef.1p65.ss[,4], FUN=function(x) x >= abs(coef.1p651.ss[1,4])))
p.coef1p65.5<-mean(sapply(coef.1p65.ss[,5], FUN=function(x) x >= abs(coef.1p651.ss[1,5])))
p.coef1p65.6<-mean(sapply(coef.1p65.ss[,6], FUN=function(x) x >= abs(coef.1p651.ss[1,6])))
p.coef1p65.7<-mean(sapply(coef.1p65.ss[,7], FUN=function(x) x >= abs(coef.1p651.ss[1,7])))
p.coef1p65.8<-mean(sapply(coef.1p65.ss[,8], FUN=function(x) x >= abs(coef.1p651.ss[1,8])))
p.coef1p65.9<-mean(sapply(coef.1p65.ss[,9], FUN=function(x) x >= abs(coef.1p651.ss[1,9])))
p.coef1p65.10<-mean(sapply(coef.1p65.ss[,10], FUN=function(x) x >= abs(coef.1p651.ss[1,10])))

pvals_rand[5,]<-c(p.coef1p65.1, p.coef1p65.2, p.coef1p65.3, p.coef1p65.4, p.coef1p65.5,
                  p.coef1p65.6, p.coef1p65.7, p.coef1p65.8,  p.coef1p65.9, p.coef1p65.10)

discarded_iters[5,]<-10000-nrow(coef.1p65.ss)


#### Cutoff 1.85
bin_df1p85_ss<-glm(data=df, victory~onep85cut+op85_par+proximity +  polity21 + 
                     capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.1p85.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$onep85cut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.1p85.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.1p85.ss<-cbind(coef.1p85.ss, fail)
coef.1p85.ss<-subset(coef.1p85.ss, coef.1p85.ss[,11]==0)

coef.1p851.ss<-matrix(NA, 1, 10)
coef.1p851.ss[1,]<-summary(bin_df1p85_ss)$coefficients[,3]

p.coef1p85.1<-mean(sapply(coef.1p85.ss[,1], FUN=function(x) x >= abs(coef.1p851.ss[1,1])))
p.coef1p85.2<-mean(sapply(coef.1p85.ss[,2], FUN=function(x) x >= abs(coef.1p851.ss[1,2])))
p.coef1p85.3<-mean(sapply(coef.1p85.ss[,3], FUN=function(x) x >= abs(coef.1p851.ss[1,3])))
p.coef1p85.4<-mean(sapply(coef.1p85.ss[,4], FUN=function(x) x >= abs(coef.1p851.ss[1,4])))
p.coef1p85.5<-mean(sapply(coef.1p85.ss[,5], FUN=function(x) x >= abs(coef.1p851.ss[1,5])))
p.coef1p85.6<-mean(sapply(coef.1p85.ss[,6], FUN=function(x) x >= abs(coef.1p851.ss[1,6])))
p.coef1p85.7<-mean(sapply(coef.1p85.ss[,7], FUN=function(x) x >= abs(coef.1p851.ss[1,7])))
p.coef1p85.8<-mean(sapply(coef.1p85.ss[,8], FUN=function(x) x >= abs(coef.1p851.ss[1,8])))
p.coef1p85.9<-mean(sapply(coef.1p85.ss[,9], FUN=function(x) x >= abs(coef.1p851.ss[1,9])))
p.coef1p85.10<-mean(sapply(coef.1p85.ss[,10], FUN=function(x) x >= abs(coef.1p851.ss[1,10])))

pvals_rand[6,]<-c(p.coef1p85.1, p.coef1p85.2, p.coef1p85.3, p.coef1p85.4, p.coef1p85.5,
                  p.coef1p85.6, p.coef1p85.7, p.coef1p85.8,  p.coef1p85.9, p.coef1p85.10)

discarded_iters[6,]<-10000-nrow(coef.1p85.ss)

#### Cutoff 2
bin_df2_ss<-glm(data=df, victory~twocut+two_par+proximity +  polity21 + 
                  capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.2.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$twocut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.2.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.2.ss<-cbind(coef.2.ss, fail)
coef.2.ss<-subset(coef.2.ss, coef.2.ss[,11]==0)

coef.21.ss<-matrix(NA, 1, 10)
coef.21.ss[1,]<-summary(bin_df2_ss)$coefficients[,3]

p.coef2.1<-mean(sapply(coef.2.ss[,1], FUN=function(x) x >= abs(coef.21.ss[1,1])))
p.coef2.2<-mean(sapply(coef.2.ss[,2], FUN=function(x) x >= abs(coef.21.ss[1,2])))
p.coef2.3<-mean(sapply(coef.2.ss[,3], FUN=function(x) x >= abs(coef.21.ss[1,3])))
p.coef2.4<-mean(sapply(coef.2.ss[,4], FUN=function(x) x >= abs(coef.21.ss[1,4])))
p.coef2.5<-mean(sapply(coef.2.ss[,5], FUN=function(x) x >= abs(coef.21.ss[1,5])))
p.coef2.6<-mean(sapply(coef.2.ss[,6], FUN=function(x) x >= abs(coef.21.ss[1,6])))
p.coef2.7<-mean(sapply(coef.2.ss[,7], FUN=function(x) x >= abs(coef.21.ss[1,7])))
p.coef2.8<-mean(sapply(coef.2.ss[,8], FUN=function(x) x >= abs(coef.21.ss[1,8])))
p.coef2.9<-mean(sapply(coef.2.ss[,9], FUN=function(x) x >= abs(coef.21.ss[1,9])))
p.coef2.10<-mean(sapply(coef.2.ss[,10], FUN=function(x) x >= abs(coef.21.ss[1,10])))

pvals_rand[7,]<-c(p.coef2.1, p.coef2.2, p.coef2.3, p.coef2.4, p.coef2.5,
                  p.coef2.6, p.coef2.7, p.coef2.8,  p.coef2.9, p.coef2.10)

discarded_iters[7,]<-10000-nrow(coef.2.ss)


####
bin_df25_ss<-glm(data=df, victory~twop5cut+tp5_par+proximity +  polity21 + 
                   capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.2p5.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$twop5cut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.2p5.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.2p5.ss<-cbind(coef.2p5.ss, fail)
coef.2p5.ss<-subset(coef.2p5.ss, coef.2p5.ss[,11]==0)

coef.2p51.ss<-matrix(NA, 1, 10)
coef.2p51.ss[1,]<-summary(bin_df25_ss)$coefficients[,3]

p.coef2p5.1<-mean(sapply(coef.2p5.ss[,1], FUN=function(x) x >= abs(coef.2p51.ss[1,1])))
p.coef2p5.2<-mean(sapply(coef.2.ss[,2], FUN=function(x) x >= abs(coef.2p51.ss[1,2])))
p.coef2p5.3<-mean(sapply(coef.2p5.ss[,3], FUN=function(x) x >= abs(coef.2p51.ss[1,3])))
p.coef2p5.4<-mean(sapply(coef.2p5.ss[,4], FUN=function(x) x >= abs(coef.2p51.ss[1,4])))
p.coef2p5.5<-mean(sapply(coef.2p5.ss[,5], FUN=function(x) x >= abs(coef.2p51.ss[1,5])))
p.coef2p5.6<-mean(sapply(coef.2p5.ss[,6], FUN=function(x) x >= abs(coef.2p51.ss[1,6])))
p.coef2p5.7<-mean(sapply(coef.2p5.ss[,7], FUN=function(x) x >= abs(coef.2p51.ss[1,7])))
p.coef2p5.8<-mean(sapply(coef.2p5.ss[,8], FUN=function(x) x >= abs(coef.2p51.ss[1,8])))
p.coef2p5.9<-mean(sapply(coef.2p5.ss[,9], FUN=function(x) x >= abs(coef.2p51.ss[1,9])))
p.coef2p5.10<-mean(sapply(coef.2p5.ss[,10], FUN=function(x) x >= abs(coef.2p51.ss[1,10])))


pvals_rand[8,]<-c(p.coef2p5.1, p.coef2p5.2, p.coef2p5.3, p.coef2p5.4, p.coef2p5.5,
                  p.coef2p5.6, p.coef2p5.7, p.coef2p5.8,  p.coef2p5.9, p.coef2p5.10)

discarded_iters[8,]<-10000-nrow(coef.2p5.ss)


#### Cutoff 3
bin_df3_ss<-glm(data=df, victory~threecut+three_par+proximity +  polity21 + 
                  capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.3.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$threecut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.3.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.3.ss<-cbind(coef.3.ss, fail)
coef.3.ss<-subset(coef.3.ss, coef.3.ss[,11]==0)


coef.31.ss<-matrix(NA, 1, 10)
coef.31.ss[1,]<-summary(bin_df3_ss)$coefficients[,3]

p.coef3.1<-mean(sapply(coef.3.ss[,1], FUN=function(x) x >= abs(coef.31.ss[1,1])))
p.coef3.2<-mean(sapply(coef.3.ss[,2], FUN=function(x) x >= abs(coef.31.ss[1,2])))
p.coef3.3<-mean(sapply(coef.3.ss[,3], FUN=function(x) x >= abs(coef.31.ss[1,3])))
p.coef3.4<-mean(sapply(coef.3.ss[,4], FUN=function(x) x >= abs(coef.31.ss[1,4])))
p.coef3.5<-mean(sapply(coef.3.ss[,5], FUN=function(x) x >= abs(coef.31.ss[1,5])))
p.coef3.6<-mean(sapply(coef.3.ss[,6], FUN=function(x) x >= abs(coef.31.ss[1,6])))
p.coef3.7<-mean(sapply(coef.3.ss[,7], FUN=function(x) x >= abs(coef.31.ss[1,7])))
p.coef3.8<-mean(sapply(coef.3.ss[,8], FUN=function(x) x >= abs(coef.31.ss[1,8])))
p.coef3.9<-mean(sapply(coef.3.ss[,9], FUN=function(x) x >= abs(coef.31.ss[1,9])))
p.coef3.10<-mean(sapply(coef.3.ss[,10], FUN=function(x) x >= abs(coef.31.ss[1,10])))

pvals_rand[9,]<-c(p.coef3.1, p.coef3.2, p.coef3.3, p.coef3.4, p.coef3.5,
                  p.coef3.6, p.coef3.7, p.coef3.8,  p.coef3.9, p.coef3.10)

discarded_iters[9,]<-10000-nrow(coef.3.ss)

#### Cutoff 3.4
bin_df3p4_ss<-glm(data=df, victory~threep4cut+ threep4_par + proximity +  polity21 + 
                    capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.3p4.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$threep4cut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.3p4.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.3p4.ss<-cbind(coef.3p4.ss, fail)
coef.3p4.ss<-subset(coef.3p4.ss, coef.3p4.ss[,11]==0)

coef.3p41.ss<-matrix(NA, 1, 10)
coef.3p41.ss[1,]<-summary(bin_df3p4_ss)$coefficients[,3]

p.coef3p4.1<-mean(sapply(coef.3p4.ss[,1], FUN=function(x) x >= abs(coef.3p41.ss[1,1])))
p.coef3p4.2<-mean(sapply(coef.3p4.ss[,2], FUN=function(x) x >= abs(coef.3p41.ss[1,2])))
p.coef3p4.3<-mean(sapply(coef.3p4.ss[,3], FUN=function(x) x >= abs(coef.3p41.ss[1,3])))
p.coef3p4.4<-mean(sapply(coef.3p4.ss[,4], FUN=function(x) x >= abs(coef.3p41.ss[1,4])))
p.coef3p4.5<-mean(sapply(coef.3p4.ss[,5], FUN=function(x) x >= abs(coef.3p41.ss[1,5])))
p.coef3p4.6<-mean(sapply(coef.3p4.ss[,6], FUN=function(x) x >= abs(coef.3p41.ss[1,6])))
p.coef3p4.7<-mean(sapply(coef.3p4.ss[,7], FUN=function(x) x >= abs(coef.3p41.ss[1,7])))
p.coef3p4.8<-mean(sapply(coef.3p4.ss[,8], FUN=function(x) x >= abs(coef.3p41.ss[1,8])))
p.coef3p4.9<-mean(sapply(coef.3p4.ss[,9], FUN=function(x) x >= abs(coef.3p41.ss[1,9])))
p.coef3p4.10<-mean(sapply(coef.3p4.ss[,10], FUN=function(x) x >= abs(coef.3p41.ss[1,10])))

pvals_rand[10,]<-c(p.coef3p4.1, p.coef3p4.2, p.coef3p4.3, p.coef3p4.4, p.coef3p4.5,
                   p.coef3p4.6, p.coef3p4.7, p.coef3p4.8,  p.coef3p4.9, p.coef3p4.10)

discarded_iters[10,]<-10000-nrow(coef.3p4.ss)


#### Cutoff 3.5
bin_df3p5_ss<-glm(data=df, victory~threep5cut+ threep5_par + proximity +  polity21 + 
                    capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.3p5.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$threep5cut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.3p5.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.3p5.ss<-cbind(coef.3p5.ss, fail)
coef.3p5.ss<-subset(coef.3p5.ss, coef.3p5.ss[,11]==0)

coef.3p51.ss<-matrix(NA, 1, 10)
coef.3p51.ss[1,]<-summary(bin_df3p5_ss)$coefficients[,3]

p.coef3p5.1<-mean(sapply(coef.3p5.ss[,1], FUN=function(x) x >= abs(coef.3p51.ss[1,1])))
p.coef3p5.2<-mean(sapply(coef.3p5.ss[,2], FUN=function(x) x >= abs(coef.3p51.ss[1,2])))
p.coef3p5.3<-mean(sapply(coef.3p5.ss[,3], FUN=function(x) x >= abs(coef.3p51.ss[1,3])))
p.coef3p5.4<-mean(sapply(coef.3p5.ss[,4], FUN=function(x) x >= abs(coef.3p51.ss[1,4])))
p.coef3p5.5<-mean(sapply(coef.3p5.ss[,5], FUN=function(x) x >= abs(coef.3p51.ss[1,5])))
p.coef3p5.6<-mean(sapply(coef.3p5.ss[,6], FUN=function(x) x >= abs(coef.3p51.ss[1,6])))
p.coef3p5.7<-mean(sapply(coef.3p5.ss[,7], FUN=function(x) x >= abs(coef.3p51.ss[1,7])))
p.coef3p5.8<-mean(sapply(coef.3p5.ss[,8], FUN=function(x) x >= abs(coef.3p51.ss[1,8])))
p.coef3p5.9<-mean(sapply(coef.3p5.ss[,9], FUN=function(x) x >= abs(coef.3p51.ss[1,9])))
p.coef3p5.10<-mean(sapply(coef.3p5.ss[,10], FUN=function(x) x >= abs(coef.3p51.ss[1,10])))

pvals_rand[11,]<-c(p.coef3p5.1, p.coef3p5.2, p.coef3p5.3, p.coef3p5.4, p.coef3p5.5,
                   p.coef3p5.6, p.coef3p5.7, p.coef3p5.8,  p.coef3p5.9, p.coef3p5.10)

discarded_iters[11,]<-10000-nrow(coef.3p5.ss)


#### 4 CUTOFF
bin_dffour_ss<-glm(data=df, victory~fourcut+four_par+proximity +  polity21 + 
                     capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.four.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)


reps<-sum(df$fourcut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.four.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.four.ss<-cbind(coef.four.ss, fail)
coef.four.ss<-subset(coef.four.ss, coef.four.ss[,11]==0)

coef.four1.ss<-matrix(NA, 1, 10)
coef.four1.ss[1,]<-summary(bin_dffour_ss)$coefficients[,3]

p.coeffour.1<-mean(sapply(coef.four.ss[,1], FUN=function(x) x >= abs(coef.four1.ss[1,1])))
p.coeffour.2<-mean(sapply(coef.four.ss[,2], FUN=function(x) x >= abs(coef.four1.ss[1,2])))
p.coeffour.3<-mean(sapply(coef.four.ss[,3], FUN=function(x) x >= abs(coef.four1.ss[1,3])))
p.coeffour.4<-mean(sapply(coef.four.ss[,4], FUN=function(x) x >= abs(coef.four1.ss[1,4])))
p.coeffour.5<-mean(sapply(coef.four.ss[,5], FUN=function(x) x >= abs(coef.four1.ss[1,5])))
p.coeffour.6<-mean(sapply(coef.four.ss[,6], FUN=function(x) x >= abs(coef.four1.ss[1,6])))
p.coeffour.7<-mean(sapply(coef.four.ss[,7], FUN=function(x) x >= abs(coef.four1.ss[1,7])))
p.coeffour.8<-mean(sapply(coef.four.ss[,8], FUN=function(x) x >= abs(coef.four1.ss[1,8])))
p.coeffour.9<-mean(sapply(coef.four.ss[,9], FUN=function(x) x >= abs(coef.four1.ss[1,9])))
p.coeffour.10<-mean(sapply(coef.four.ss[,10], FUN=function(x) x >= abs(coef.four1.ss[1,10])))

pvals_rand[12,]<-c(p.coeffour.1, p.coeffour.2, p.coeffour.3, p.coeffour.4, p.coeffour.5,
                   p.coeffour.6, p.coeffour.7, p.coeffour.8,  p.coeffour.9, p.coeffour.10)

discarded_iters[12,]<-10000-nrow(coef.four.ss)

##### 7 CUTOFF
bin_df7_ss<-glm(data=df, victory~sevencut+seven_par+proximity +  polity21 + 
                  capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.7.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$sevencut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.7.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.7.ss<-cbind(coef.7.ss, fail)
coef.7.ss<-subset(coef.7.ss, coef.7.ss[,11]==0)

coef.71.ss<-matrix(NA, 1, 10)
coef.71.ss[1,]<-summary(bin_df7_ss)$coefficients[,3]

p.coef7.1<-mean(sapply(coef.7.ss[,1], FUN=function(x) x >= abs(coef.71.ss[1,1])))
p.coef7.2<-mean(sapply(coef.7.ss[,2], FUN=function(x) x >= abs(coef.71.ss[1,2])))
p.coef7.3<-mean(sapply(coef.7.ss[,3], FUN=function(x) x >= abs(coef.71.ss[1,3])))
p.coef7.4<-mean(sapply(coef.7.ss[,4], FUN=function(x) x >= abs(coef.71.ss[1,4])))
p.coef7.5<-mean(sapply(coef.7.ss[,5], FUN=function(x) x >= abs(coef.71.ss[1,5])))
p.coef7.6<-mean(sapply(coef.7.ss[,6], FUN=function(x) x >= abs(coef.71.ss[1,6])))
p.coef7.7<-mean(sapply(coef.7.ss[,7], FUN=function(x) x >= abs(coef.71.ss[1,7])))
p.coef7.8<-mean(sapply(coef.7.ss[,8], FUN=function(x) x >= abs(coef.71.ss[1,8])))
p.coef7.9<-mean(sapply(coef.7.ss[,9], FUN=function(x) x >= abs(coef.71.ss[1,9])))
p.coef7.10<-mean(sapply(coef.7.ss[,10], FUN=function(x) x >= abs(coef.71.ss[1,10])))

pvals_rand[13,]<-c(p.coef7.1, p.coef7.2, p.coef7.3, p.coef7.4, p.coef7.5,
                   p.coef7.6, p.coef7.7, p.coef7.8,  p.coef7.9, p.coef7.10)

discarded_iters[13,]<-10000-nrow(coef.7.ss)


##### 9 CUTOFF
bin_df9_ss<-glm(data=df, victory~ninecut+nine_par+proximity +  polity21 + 
                  capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.9.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$ninecut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.9.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.9.ss<-cbind(coef.9.ss, fail)
coef.9.ss<-subset(coef.9.ss, coef.9.ss[,11]==0)

coef.91.ss<-matrix(NA, 1, 10)
coef.91.ss[1,]<-summary(bin_df9_ss)$coefficients[,3]

p.coef9.1<-mean(sapply(coef.9.ss[,1], FUN=function(x) x >= abs(coef.91.ss[1,1])))
p.coef9.2<-mean(sapply(coef.9.ss[,2], FUN=function(x) x >= abs(coef.91.ss[1,2])))
p.coef9.3<-mean(sapply(coef.9.ss[,3], FUN=function(x) x >= abs(coef.91.ss[1,3])))
p.coef9.4<-mean(sapply(coef.9.ss[,4], FUN=function(x) x >= abs(coef.91.ss[1,4])))
p.coef9.5<-mean(sapply(coef.9.ss[,5], FUN=function(x) x >= abs(coef.91.ss[1,5])))
p.coef9.6<-mean(sapply(coef.9.ss[,6], FUN=function(x) x >= abs(coef.91.ss[1,6])))
p.coef9.7<-mean(sapply(coef.9.ss[,7], FUN=function(x) x >= abs(coef.91.ss[1,7])))
p.coef9.8<-mean(sapply(coef.9.ss[,8], FUN=function(x) x >= abs(coef.91.ss[1,8])))
p.coef9.9<-mean(sapply(coef.9.ss[,9], FUN=function(x) x >= abs(coef.91.ss[1,9])))
p.coef9.10<-mean(sapply(coef.9.ss[,10], FUN=function(x) x >= abs(coef.91.ss[1,10])))

pvals_rand[14,]<-c(p.coef9.1, p.coef9.2, p.coef9.3, p.coef9.4, p.coef9.5,
                   p.coef9.6, p.coef9.7, p.coef9.8,  p.coef9.9, p.coef9.10)

discarded_iters[14,]<-10000-nrow(coef.9.ss)


### CUTOFF 10 
bin_df10_ss<-glm(data=df, victory~tencut+ten_par+proximity +  polity21 + 
                   capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.10.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$tencut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.10.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.10.ss<-cbind(coef.10.ss, fail)
coef.10.ss<-subset(coef.10.ss, coef.10.ss[,11]==0)

coef.101.ss<-matrix(NA, 1, 10)
coef.101.ss[1,]<-summary(bin_df10_ss)$coefficients[,3]

p.coef10.1<-mean(sapply(coef.10.ss[,1], FUN=function(x) x >= abs(coef.101.ss[1,1])))
p.coef10.2<-mean(sapply(coef.10.ss[,2], FUN=function(x) x >= abs(coef.101.ss[1,2])))
p.coef10.3<-mean(sapply(coef.10.ss[,3], FUN=function(x) x >= abs(coef.101.ss[1,3])))
p.coef10.4<-mean(sapply(coef.10.ss[,4], FUN=function(x) x >= abs(coef.101.ss[1,4])))
p.coef10.5<-mean(sapply(coef.10.ss[,5], FUN=function(x) x >= abs(coef.101.ss[1,5])))
p.coef10.6<-mean(sapply(coef.10.ss[,6], FUN=function(x) x >= abs(coef.101.ss[1,6])))
p.coef10.7<-mean(sapply(coef.10.ss[,7], FUN=function(x) x >= abs(coef.101.ss[1,7])))
p.coef10.8<-mean(sapply(coef.10.ss[,8], FUN=function(x) x >= abs(coef.101.ss[1,8])))
p.coef10.9<-mean(sapply(coef.10.ss[,9], FUN=function(x) x >= abs(coef.101.ss[1,9])))
p.coef10.10<-mean(sapply(coef.10.ss[,10], FUN=function(x) x >= abs(coef.101.ss[1,10])))

pvals_rand[15,]<-c(p.coef10.1, p.coef10.2, p.coef10.3, p.coef10.4, p.coef10.5,
                   p.coef10.6, p.coef10.7, p.coef10.8,  p.coef10.9, p.coef10.10)

discarded_iters[15,]<-10000-nrow(coef.10.ss)


### CUTOFF 11 
bin_df11_ss<-glm(data=df, victory~elevencut+elevencut_par+proximity +  polity21 + 
                   capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.11.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$elevencut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.11.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.11.ss<-cbind(coef.11.ss, fail)
coef.11.ss<-subset(coef.11.ss, coef.11.ss[,11]==0)

coef.111.ss<-matrix(NA, 1, 10)
coef.111.ss[1,]<-summary(bin_df11_ss)$coefficients[,3]

p.coef11.1<-mean(sapply(coef.11.ss[,1], FUN=function(x) x >= abs(coef.111.ss[1,1])))
p.coef11.2<-mean(sapply(coef.11.ss[,2], FUN=function(x) x >= abs(coef.111.ss[1,2])))
p.coef11.3<-mean(sapply(coef.11.ss[,3], FUN=function(x) x >= abs(coef.111.ss[1,3])))
p.coef11.4<-mean(sapply(coef.11.ss[,4], FUN=function(x) x >= abs(coef.111.ss[1,4])))
p.coef11.5<-mean(sapply(coef.11.ss[,5], FUN=function(x) x >= abs(coef.111.ss[1,5])))
p.coef11.6<-mean(sapply(coef.11.ss[,6], FUN=function(x) x >= abs(coef.111.ss[1,6])))
p.coef11.7<-mean(sapply(coef.11.ss[,7], FUN=function(x) x >= abs(coef.111.ss[1,7])))
p.coef11.8<-mean(sapply(coef.11.ss[,8], FUN=function(x) x >= abs(coef.111.ss[1,8])))
p.coef11.9<-mean(sapply(coef.11.ss[,9], FUN=function(x) x >= abs(coef.111.ss[1,9])))
p.coef11.10<-mean(sapply(coef.11.ss[,10], FUN=function(x) x >= abs(coef.111.ss[1,10])))

pvals_rand[16,]<-c(p.coef11.1, p.coef11.2, p.coef11.3, p.coef11.4, p.coef11.5,
                   p.coef11.6, p.coef11.7, p.coef11.8,  p.coef11.9, p.coef11.10)

discarded_iters[16,]<-10000-nrow(coef.11.ss)

### CUTOFF 20
bin_df20_ss<-glm(data=df, victory~twentycut+twenty_par+proximity +  polity21 + 
                   capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.20.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$twentycut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.20.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.20.ss<-cbind(coef.20.ss, fail)
coef.20.ss<-subset(coef.20.ss, coef.20.ss[,11]==0)

coef.201.ss<-matrix(NA, 1, 10)
coef.201.ss[1,]<-summary(bin_df20_ss)$coefficients[,3]

p.coef20.1<-mean(sapply(coef.20.ss[,1], FUN=function(x) x >= abs(coef.201.ss[1,1])))
p.coef20.2<-mean(sapply(coef.20.ss[,2], FUN=function(x) x >= abs(coef.201.ss[1,2])))
p.coef20.3<-mean(sapply(coef.20.ss[,3], FUN=function(x) x >= abs(coef.201.ss[1,3])))
p.coef20.4<-mean(sapply(coef.20.ss[,4], FUN=function(x) x >= abs(coef.201.ss[1,4])))
p.coef20.5<-mean(sapply(coef.20.ss[,5], FUN=function(x) x >= abs(coef.201.ss[1,5])))
p.coef20.6<-mean(sapply(coef.20.ss[,6], FUN=function(x) x >= abs(coef.201.ss[1,6])))
p.coef20.7<-mean(sapply(coef.20.ss[,7], FUN=function(x) x >= abs(coef.201.ss[1,7])))
p.coef20.8<-mean(sapply(coef.20.ss[,8], FUN=function(x) x >= abs(coef.201.ss[1,8])))
p.coef20.9<-mean(sapply(coef.20.ss[,9], FUN=function(x) x >= abs(coef.201.ss[1,9])))
p.coef20.10<-mean(sapply(coef.20.ss[,10], FUN=function(x) x >= abs(coef.201.ss[1,10])))

pvals_rand[17,]<-c(p.coef20.1, p.coef20.2, p.coef20.3, p.coef20.4, p.coef20.5,
                   p.coef20.6, p.coef20.7, p.coef20.8,  p.coef20.9, p.coef20.10)

discarded_iters[17,]<-10000-nrow(coef.20.ss)

### Cutoff 28
bin_df28_ss<-glm(data=df, victory~twentyeightcut+twentyeight_par+proximity +  polity21 + 
                   capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.28.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$twentyeightcut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.28.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.28.ss<-cbind(coef.28.ss, fail)
coef.28.ss<-subset(coef.28.ss, coef.28.ss[,11]==0)

coef.281.ss<-matrix(NA, 1, 10)
coef.281.ss[1,]<-summary(bin_df28_ss)$coefficients[,3]

p.coef28.1<-mean(sapply(coef.28.ss[,1], FUN=function(x) x >= abs(coef.281.ss[1,1])))
p.coef28.2<-mean(sapply(coef.28.ss[,2], FUN=function(x) x >= abs(coef.281.ss[1,2])))
p.coef28.3<-mean(sapply(coef.28.ss[,3], FUN=function(x) x >= abs(coef.281.ss[1,3])))
p.coef28.4<-mean(sapply(coef.28.ss[,4], FUN=function(x) x >= abs(coef.281.ss[1,4])))
p.coef28.5<-mean(sapply(coef.28.ss[,5], FUN=function(x) x >= abs(coef.281.ss[1,5])))
p.coef28.6<-mean(sapply(coef.28.ss[,6], FUN=function(x) x >= abs(coef.281.ss[1,6])))
p.coef28.7<-mean(sapply(coef.28.ss[,7], FUN=function(x) x >= abs(coef.281.ss[1,7])))
p.coef28.8<-mean(sapply(coef.28.ss[,8], FUN=function(x) x >= abs(coef.281.ss[1,8])))
p.coef28.9<-mean(sapply(coef.28.ss[,9], FUN=function(x) x >= abs(coef.281.ss[1,9])))
p.coef28.10<-mean(sapply(coef.28.ss[,10], FUN=function(x) x >= abs(coef.281.ss[1,10])))

pvals_rand[18,]<-c(p.coef28.1, p.coef28.2, p.coef28.3, p.coef28.4, p.coef28.5,
                   p.coef28.6, p.coef28.7, p.coef28.8,  p.coef28.9, p.coef28.10)

discarded_iters[18,]<-10000-nrow(coef.28.ss)

### Cutoff 30
bin_df30_ss<-glm(data=df, victory~thirtycut+thirty_par+proximity +  polity21 + 
                   capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.30.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$thirtycut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.30.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.30.ss<-cbind(coef.30.ss, fail)
coef.30.ss<-subset(coef.30.ss, coef.30.ss[,11]==0)

coef.301.ss<-matrix(NA, 1, 10)
coef.301.ss[1,]<-summary(bin_df30_ss)$coefficients[,3]

p.coef30.1<-mean(sapply(coef.30.ss[,1], FUN=function(x) x >= abs(coef.301.ss[1,1])))
p.coef30.2<-mean(sapply(coef.30.ss[,2], FUN=function(x) x >= abs(coef.301.ss[1,2])))
p.coef30.3<-mean(sapply(coef.30.ss[,3], FUN=function(x) x >= abs(coef.301.ss[1,3])))
p.coef30.4<-mean(sapply(coef.30.ss[,4], FUN=function(x) x >= abs(coef.301.ss[1,4])))
p.coef30.5<-mean(sapply(coef.30.ss[,5], FUN=function(x) x >= abs(coef.301.ss[1,5])))
p.coef30.6<-mean(sapply(coef.30.ss[,6], FUN=function(x) x >= abs(coef.301.ss[1,6])))
p.coef30.7<-mean(sapply(coef.30.ss[,7], FUN=function(x) x >= abs(coef.301.ss[1,7])))
p.coef30.8<-mean(sapply(coef.30.ss[,8], FUN=function(x) x >= abs(coef.301.ss[1,8])))
p.coef30.9<-mean(sapply(coef.30.ss[,9], FUN=function(x) x >= abs(coef.301.ss[1,9])))
p.coef30.10<-mean(sapply(coef.30.ss[,10], FUN=function(x) x >= abs(coef.301.ss[1,10])))

pvals_rand[19,]<-c(p.coef30.1, p.coef30.2, p.coef30.3, p.coef30.4, p.coef30.5,
                   p.coef30.6, p.coef30.7, p.coef30.8,  p.coef30.9, p.coef30.10)

discarded_iters[19,]<-10000-nrow(coef.30.ss)

#### CUTOFF 40
bin_df40_ss<-glm(data=df, victory~fortycut+forty_par+proximity +  polity21 + 
                   capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.40.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$fortycut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.40.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.40.ss<-cbind(coef.40.ss, fail)
coef.40.ss<-subset(coef.40.ss, coef.40.ss[,11]==0)

coef.401.ss<-matrix(NA, 1, 10)
coef.401.ss[1,]<-summary(bin_df40_ss)$coefficients[,3]

p.coef40.1<-mean(sapply(coef.40.ss[,1], FUN=function(x) x >= abs(coef.401.ss[1,1])))
p.coef40.2<-mean(sapply(coef.40.ss[,2], FUN=function(x) x >= abs(coef.401.ss[1,2])))
p.coef40.3<-mean(sapply(coef.40.ss[,3], FUN=function(x) x >= abs(coef.401.ss[1,3])))
p.coef40.4<-mean(sapply(coef.40.ss[,4], FUN=function(x) x >= abs(coef.401.ss[1,4])))
p.coef40.5<-mean(sapply(coef.40.ss[,5], FUN=function(x) x >= abs(coef.401.ss[1,5])))
p.coef40.6<-mean(sapply(coef.40.ss[,6], FUN=function(x) x >= abs(coef.401.ss[1,6])))
p.coef40.7<-mean(sapply(coef.40.ss[,7], FUN=function(x) x >= abs(coef.401.ss[1,7])))
p.coef40.8<-mean(sapply(coef.40.ss[,8], FUN=function(x) x >= abs(coef.401.ss[1,8])))
p.coef40.9<-mean(sapply(coef.40.ss[,9], FUN=function(x) x >= abs(coef.401.ss[1,9])))
p.coef40.10<-mean(sapply(coef.40.ss[,10], FUN=function(x) x >= abs(coef.401.ss[1,10])))

pvals_rand[20,]<-c(p.coef40.1, p.coef40.2, p.coef40.3, p.coef40.4, p.coef40.5,
                   p.coef40.6, p.coef40.7, p.coef40.8,  p.coef40.9, p.coef40.10)

discarded_iters[20,]<-10000-nrow(coef.40.ss)

#### CUTOFF 50
bin_df50_ss<-glm(data=df, victory~fiftycut+fifty_par+proximity +  polity21 + 
                   capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.50.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$fiftycut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.50.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.50.ss<-cbind(coef.50.ss, fail)
coef.50.ss<-subset(coef.50.ss, coef.50.ss[,11]==0)

coef.501.ss<-matrix(NA, 1, 10)
coef.501.ss[1,]<-summary(bin_df50_ss)$coefficients[,3]

p.coef50.1<-mean(sapply(coef.50.ss[,1], FUN=function(x) x >= abs(coef.501.ss[1,1])))
p.coef50.2<-mean(sapply(coef.50.ss[,2], FUN=function(x) x >= abs(coef.501.ss[1,2])))
p.coef50.3<-mean(sapply(coef.50.ss[,3], FUN=function(x) x >= abs(coef.501.ss[1,3])))
p.coef50.4<-mean(sapply(coef.50.ss[,4], FUN=function(x) x >= abs(coef.501.ss[1,4])))
p.coef50.5<-mean(sapply(coef.50.ss[,5], FUN=function(x) x >= abs(coef.501.ss[1,5])))
p.coef50.6<-mean(sapply(coef.50.ss[,6], FUN=function(x) x >= abs(coef.501.ss[1,6])))
p.coef50.7<-mean(sapply(coef.50.ss[,7], FUN=function(x) x >= abs(coef.501.ss[1,7])))
p.coef50.8<-mean(sapply(coef.50.ss[,8], FUN=function(x) x >= abs(coef.501.ss[1,8])))
p.coef50.9<-mean(sapply(coef.50.ss[,9], FUN=function(x) x >= abs(coef.501.ss[1,9])))
p.coef50.10<-mean(sapply(coef.50.ss[,10], FUN=function(x) x >= abs(coef.501.ss[1,10])))

pvals_rand[21,]<-c(p.coef50.1, p.coef50.2, p.coef50.3, p.coef50.4, p.coef50.5,
                   p.coef50.6, p.coef50.7, p.coef50.8,  p.coef50.9, p.coef50.10)

discarded_iters[21,]<-10000-nrow(coef.50.ss)


### CUTOFF 100
bin_df100_ss<-glm(data=df, victory~hundredcut+hundred_par+proximity +  polity21 + 
                    capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.100.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$hundredcut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.100.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.100.ss<-cbind(coef.100.ss, fail)
coef.100.ss<-subset(coef.100.ss, coef.100.ss[,11]==0)

coef.1001.ss<-matrix(NA, 1, 10)
coef.1001.ss[1,]<-summary(bin_df100_ss)$coefficients[,3]

p.coef100.1<-mean(sapply(coef.100.ss[,1], FUN=function(x) x >= abs(coef.1001.ss[1,1])))
p.coef100.2<-mean(sapply(coef.100.ss[,2], FUN=function(x) x >= abs(coef.1001.ss[1,2])))
p.coef100.3<-mean(sapply(coef.100.ss[,3], FUN=function(x) x >= abs(coef.1001.ss[1,3])))
p.coef100.4<-mean(sapply(coef.100.ss[,4], FUN=function(x) x >= abs(coef.1001.ss[1,4])))
p.coef100.5<-mean(sapply(coef.100.ss[,5], FUN=function(x) x >= abs(coef.1001.ss[1,5])))
p.coef100.6<-mean(sapply(coef.100.ss[,6], FUN=function(x) x >= abs(coef.1001.ss[1,6])))
p.coef100.7<-mean(sapply(coef.100.ss[,7], FUN=function(x) x >= abs(coef.1001.ss[1,7])))
p.coef100.8<-mean(sapply(coef.100.ss[,8], FUN=function(x) x >= abs(coef.1001.ss[1,8])))
p.coef100.9<-mean(sapply(coef.100.ss[,9], FUN=function(x) x >= abs(coef.1001.ss[1,9])))
p.coef100.10<-mean(sapply(coef.100.ss[,10], FUN=function(x) x >= abs(coef.1001.ss[1,10])))

pvals_rand[22,]<-c(p.coef100.1, p.coef100.2, p.coef100.3, p.coef100.4, p.coef100.5,
                   p.coef100.6, p.coef100.7, p.coef100.8,  p.coef100.9, p.coef100.10)

discarded_iters[22,]<-10000-nrow(coef.100.ss)

### CUTOFF 500
bin_df500_ss<-glm(data=df, victory~fivhuncut+fivhun_par+proximity +  polity21 + 
                    capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.500.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$fivhuncut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.500.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.500.ss<-cbind(coef.500.ss, fail)
coef.500.ss<-subset(coef.500.ss, coef.500.ss[,11]==0)

coef.5001.ss<-matrix(NA, 1, 10)
coef.5001.ss[1,]<-summary(bin_df500_ss)$coefficients[,3]

p.coef500.1<-mean(sapply(coef.500.ss[,1], FUN=function(x) x >= abs(coef.5001.ss[1,1])))
p.coef500.2<-mean(sapply(coef.500.ss[,2], FUN=function(x) x >= abs(coef.5001.ss[1,2])))
p.coef500.3<-mean(sapply(coef.500.ss[,3], FUN=function(x) x >= abs(coef.5001.ss[1,3])))
p.coef500.4<-mean(sapply(coef.500.ss[,4], FUN=function(x) x >= abs(coef.5001.ss[1,4])))
p.coef500.5<-mean(sapply(coef.500.ss[,5], FUN=function(x) x >= abs(coef.5001.ss[1,5])))
p.coef500.6<-mean(sapply(coef.500.ss[,6], FUN=function(x) x >= abs(coef.5001.ss[1,6])))
p.coef500.7<-mean(sapply(coef.500.ss[,7], FUN=function(x) x >= abs(coef.5001.ss[1,7])))
p.coef500.8<-mean(sapply(coef.500.ss[,8], FUN=function(x) x >= abs(coef.5001.ss[1,8])))
p.coef500.9<-mean(sapply(coef.500.ss[,9], FUN=function(x) x >= abs(coef.5001.ss[1,9])))
p.coef500.10<-mean(sapply(coef.500.ss[,10], FUN=function(x) x >= abs(coef.5001.ss[1,10])))

pvals_rand[23,]<-c(p.coef500.1, p.coef500.2, p.coef500.3, p.coef500.4, p.coef500.5,
                   p.coef500.6, p.coef500.7, p.coef500.8,  p.coef500.9, p.coef500.10)

discarded_iters[23,]<-10000-nrow(coef.500.ss)


### CUTOFF 700
bin_df700_ss<-glm(data=df, victory~sevenhuncut+sevenhuncut_par+proximity +  polity21 + 
                    capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.700.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$sevenhuncut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.700.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.700.ss<-cbind(coef.700.ss, fail)
coef.700.ss<-subset(coef.700.ss, coef.700.ss[,11]==0)

coef.7001.ss<-matrix(NA, 1, 10)
coef.7001.ss[1,]<-summary(bin_df700_ss)$coefficients[,3]

p.coef700.1<-mean(sapply(coef.700.ss[,1], FUN=function(x) x >= abs(coef.7001.ss[1,1])))
p.coef700.2<-mean(sapply(coef.700.ss[,2], FUN=function(x) x >= abs(coef.7001.ss[1,2])))
p.coef700.3<-mean(sapply(coef.700.ss[,3], FUN=function(x) x >= abs(coef.7001.ss[1,3])))
p.coef700.4<-mean(sapply(coef.700.ss[,4], FUN=function(x) x >= abs(coef.7001.ss[1,4])))
p.coef700.5<-mean(sapply(coef.700.ss[,5], FUN=function(x) x >= abs(coef.7001.ss[1,5])))
p.coef700.6<-mean(sapply(coef.700.ss[,6], FUN=function(x) x >= abs(coef.7001.ss[1,6])))
p.coef700.7<-mean(sapply(coef.700.ss[,7], FUN=function(x) x >= abs(coef.7001.ss[1,7])))
p.coef700.8<-mean(sapply(coef.700.ss[,8], FUN=function(x) x >= abs(coef.7001.ss[1,8])))
p.coef700.9<-mean(sapply(coef.700.ss[,9], FUN=function(x) x >= abs(coef.7001.ss[1,9])))
p.coef700.10<-mean(sapply(coef.700.ss[,10], FUN=function(x) x >= abs(coef.7001.ss[1,10])))

pvals_rand[24,]<-c(p.coef700.1, p.coef700.2, p.coef700.3, p.coef700.4, p.coef700.5,
                   p.coef700.6, p.coef700.7, p.coef700.8,  p.coef700.9, p.coef700.10)

discarded_iters[24,]<-10000-nrow(coef.700.ss)


### CUTOFF 1000
bin_df1000_ss<-glm(data=df, victory~thoucut+thou_par+proximity +  polity21 + 
                     capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.1000.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$thoucut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.1000.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.1000.ss<-cbind(coef.1000.ss, fail)
coef.1000.ss<-subset(coef.1000.ss, coef.1000.ss[,11]==0)

coef.10001.ss<-matrix(NA, 1, 10)
coef.10001.ss[1,]<-summary(bin_df1000_ss)$coefficients[,3]

p.coef1000.1<-mean(sapply(coef.1000.ss[,1], FUN=function(x) x >= abs(coef.10001.ss[1,1])))
p.coef1000.2<-mean(sapply(coef.1000.ss[,2], FUN=function(x) x >= abs(coef.10001.ss[1,2])))
p.coef1000.3<-mean(sapply(coef.1000.ss[,3], FUN=function(x) x >= abs(coef.10001.ss[1,3])))
p.coef1000.4<-mean(sapply(coef.1000.ss[,4], FUN=function(x) x >= abs(coef.10001.ss[1,4])))
p.coef1000.5<-mean(sapply(coef.1000.ss[,5], FUN=function(x) x >= abs(coef.10001.ss[1,5])))
p.coef1000.6<-mean(sapply(coef.1000.ss[,6], FUN=function(x) x >= abs(coef.10001.ss[1,6])))
p.coef1000.7<-mean(sapply(coef.1000.ss[,7], FUN=function(x) x >= abs(coef.10001.ss[1,7])))
p.coef1000.8<-mean(sapply(coef.1000.ss[,8], FUN=function(x) x >= abs(coef.10001.ss[1,8])))
p.coef1000.9<-mean(sapply(coef.1000.ss[,9], FUN=function(x) x >= abs(coef.10001.ss[1,9])))
p.coef1000.10<-mean(sapply(coef.1000.ss[,10], FUN=function(x) x >= abs(coef.10001.ss[1,10])))

pvals_rand[25,]<-c(p.coef1000.1, p.coef1000.2, p.coef1000.3, p.coef1000.4, p.coef1000.5,
                   p.coef1000.6, p.coef1000.7, p.coef1000.8,  p.coef1000.9, p.coef1000.10)


discarded_iters[25,]<-10000-nrow(coef.1000.ss)


### CUTOFF 1300
bin_df1300_ss<-glm(data=df, victory~thirhuncut + thirhun_par+proximity +  polity21 + 
                     capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.1300.ss<-matrix(NA, sims, 10)
coefs<-matrix(NA, sims, 10)
set.seed(1e-08)

reps<-sum(df$thirhuncut) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has superiority
    sidea<-c(rep(0,28-reps), rep(1, reps)) 
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 28, 2)
    c[,ran]<-sample(sidea, 28, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:56,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$newcut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(newcut))
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(inf_new= ifelse(sum_ord==1, 1, 0))
    df.new$inf_new<-ifelse(df.new$newcut==1, 0, df.new$inf_new)
    df.new$par_new<-ifelse(df.new$newcut==0 & df.new$inf_new==0, 1, 0)
    
    model<-glm(data=df.new, victory~newcut+par_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.1300.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.1300.ss<-cbind(coef.1300.ss, fail)
coef.1300.ss<-subset(coef.1300.ss, coef.1300.ss[,11]==0)

coef.13001.ss<-matrix(NA, 1, 10)
coef.13001.ss[1,]<-summary(bin_df1000_ss)$coefficients[,3]

p.coef1300.1<-mean(sapply(coef.1300.ss[,1], FUN=function(x) x >= abs(coef.13001.ss[1,1])))
p.coef1300.2<-mean(sapply(coef.1300.ss[,2], FUN=function(x) x >= abs(coef.13001.ss[1,2])))
p.coef1300.3<-mean(sapply(coef.1300.ss[,3], FUN=function(x) x >= abs(coef.13001.ss[1,3])))
p.coef1300.4<-mean(sapply(coef.1300.ss[,4], FUN=function(x) x >= abs(coef.13001.ss[1,4])))
p.coef1300.5<-mean(sapply(coef.1300.ss[,5], FUN=function(x) x >= abs(coef.13001.ss[1,5])))
p.coef1300.6<-mean(sapply(coef.1300.ss[,6], FUN=function(x) x >= abs(coef.13001.ss[1,6])))
p.coef1300.7<-mean(sapply(coef.1300.ss[,7], FUN=function(x) x >= abs(coef.13001.ss[1,7])))
p.coef1300.8<-mean(sapply(coef.1300.ss[,8], FUN=function(x) x >= abs(coef.13001.ss[1,8])))
p.coef1300.9<-mean(sapply(coef.1300.ss[,9], FUN=function(x) x >= abs(coef.13001.ss[1,9])))
p.coef1300.10<-mean(sapply(coef.1300.ss[,10], FUN=function(x) x >= abs(coef.13001.ss[1,10])))

pvals_rand[26,]<-c(p.coef1300.1, p.coef1300.2, p.coef1300.3, p.coef1300.4, p.coef1300.5,
                   p.coef1300.6, p.coef1300.7, p.coef1300.8,  p.coef1300.9, p.coef1300.10)


discarded_iters[26,]<-10000-nrow(coef.1300.ss)

#### CUTOFF 1
bin_df1_ss<-glm(data=df, victory~superiority+proximity +  polity21 + 
                  capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))

sims<-10000

coef.1.ss<-matrix(NA, sims, 9)
coefs<-matrix(NA, sims, 9)
set.seed(1e-08)

reps<-sum(df$superiority) 
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:58,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sup_new = sample(0:1, 2, replace=F))
    model<-glm(data=df.new, victory~sup_new+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.1.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}


fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.1.ss<-cbind(coef.1.ss, fail)
coef.1.ss<-subset(coef.1.ss, coef.1.ss[,10]==0)

coef.11.ss<-matrix(NA, 1, 9)
coef.11.ss[1,]<-summary(bin_df1_ss)$coefficients[,3]

p.coef1.1<-mean(sapply(coef.1.ss[,1], FUN=function(x) x >= abs(coef.11.ss[1,1])))
p.coef1.2<-mean(sapply(coef.1.ss[,2], FUN=function(x) x >= abs(coef.11.ss[1,2])))
p.coef1.3<-mean(sapply(coef.1.ss[,3], FUN=function(x) x >= abs(coef.11.ss[1,3])))
p.coef1.4<-mean(sapply(coef.1.ss[,4], FUN=function(x) x >= abs(coef.11.ss[1,4])))
p.coef1.5<-mean(sapply(coef.1.ss[,5], FUN=function(x) x >= abs(coef.11.ss[1,5])))
p.coef1.6<-mean(sapply(coef.1.ss[,6], FUN=function(x) x >= abs(coef.11.ss[1,6])))
p.coef1.7<-mean(sapply(coef.1.ss[,7], FUN=function(x) x >= abs(coef.11.ss[1,7])))
p.coef1.8<-mean(sapply(coef.1.ss[,8], FUN=function(x) x >= abs(coef.11.ss[1,8])))
p.coef1.9<-mean(sapply(coef.1.ss[,9], FUN=function(x) x >= abs(coef.11.ss[1,9])))

discarded_iters_sup<-10000-nrow(coef.1.ss)
discarded_iters2<-rbind(discarded_iters_sup, discarded_iters)

pvals_sup_rand<-c(p.coef1.1, p.coef1.2, "NA", p.coef1.3, p.coef1.4, p.coef1.5,
                  p.coef1.6, p.coef1.7, p.coef1.8,  p.coef1.9)

###
PVALS<-rbind(pvals_sup_rand, pvals_rand)
colnames(PVALS)<-c("Intercept", "Superiority", "Parity", "Proximity", "Regime",
                   "Capabilities", "2nd Strike", "Population", "Violence", "Security")
rownames(PVALS)<-c("1", "1.15", "1.2",
                   "1.45", "1.55", "1.65", "1.85",
                   "2", "2.5", "3", "3.4", "3.5",
                   "4", "7", 
                   "9", "10", "11", "20", "28", "30", 
                   "40",   "50",  
                   "100", "500", "700", "1000", "1300")

## Creating a plot to visualize which thresholds have p-values below 0.1 

par(mar=c(8.1, 4.1, 4.1, 4.1), xpd=T)
plot(PVALS[,2], xaxt="n", xlab=c("Nuclear Ratio: Superiority Threshold"), 
     ylab=c("P-Value"), main=c("P-Value for Effect of Superiority, via Randomization Inference"))
axis(1, at=1:27, labels=c("1", "1.15", "1.2",
                          "1.45", "1.55", "1.65", "1.85",
                          "2", "2.5", "3", "3.4", "3.5",
                          "4", "7", 
                          "9", "10", "11", "20", "28", "30", 
                          "40",   "50",  
                          "100", "500", "700", "1000", "1300"), las=2)
abline(h=0.1, xpd=F)

PVALS<-cbind(PVALS, seq(1,27))
PVALS[,c(2, 11)]



## Checking to see how many iterations were discarded from each randomization run, to make sure not too many
par(mfrow=c(1,2))
plot(discarded_iters2, xaxt="n", xlab=c("Nuclear Ratio: Superiority Threshold"), 
     ylab=c("Iterations Discarded"), main=c("Iterations Discarded Due to Perfect Separation"), ylim=c(0,10000))
axis(1, at=1:27, labels=c("1", "1.15", "1.2",
                          "1.45", "1.55", "1.65", "1.85",
                          "2", "2.5", "3", "3.4", "3.5",
                          "4", "7", 
                          "9", "10", "11", "20", "28", "30", 
                          "40",   "50",  
                          "100", "500", "700", "1000", "1300"), las=2)
plot(discarded_iters2, xaxt="n", xlab=c("Nuclear Ratio: Superiority Threshold"), 
     ylab=c("Iterations Discarded"), main=c("Zoomed In"))
axis(1, at=1:27, labels=c("1", "1.15", "1.2",
                          "1.45", "1.55", "1.65", "1.85",
                          "2", "2.5", "3", "3.4", "3.5",
                          "4", "7", 
                          "9", "10", "11", "20", "28", "30", 
                          "40",   "50",  
                          "100", "500", "700", "1000", "1300"), las=2)


### FIGURE 5 from Appendix C - With superiority from randomization inference
dev.off()
plot(first_diffs[,7], ylim=c(-1,1), xaxt="n", xlab=c("Nuclear Ratio: Asymmetric Superiority Threshold"), 
     ylab=c("First Difference"), main=c("Difference in Prob. of Victory: Superior vs. Inferior"))
axis(1, at=1:27, labels=c("1", "1.15", "1.2",
                          "1.45", "1.55", "1.65", "1.85",
                          "2", "2.5", "3", "3.4", "3.5",
                          "4", "7", 
                          "9", "10", "11", "20", "28", "30", 
                          "40",   "50",  
                          "100", "500", "700", "1000", "1300"), las=2)
points(first_diffs[1,8], first_diffs[1,7],pch=19)
points(first_diffs[2,8], first_diffs[2,7],pch=19)
points(first_diffs[4,8], first_diffs[4,7],pch=19)
points(first_diffs[20,8], first_diffs[20,7],pch=19)
points(first_diffs[21,8], first_diffs[21,7],pch=19)
legend("bottomright", inset=c(0, -0.25),legend=c("Significant Effect, Rand-Inf"), pch=19, 
       cex=0.5, pt.cex=1)
abline(h=0, xpd=F)


#####
